//Nathan Popham
//4003-573-01: Procedural Shading
//Lab 2
//
//Strauss lighting model shader 
//
//

float sqr (float x) { return x * x; }


/*
 *  Fresnel equation
 */
float F (float angle)
{
    float kf = 1.12;
    float x = angle / (PI/2.0);
    return            
    ((1.0 / sqr(x - kf)) - (1.0 / sqr(kf))) /
    ((1.0 / sqr(1.0 - kf)) - (1.0 / sqr(kf)));
}


/*
 *  Whatever the G function is
 */
float G (float angle)
{
    float kg = 1.01;
    float x = angle / (PI/2.0);

    return
    ((1.0 / sqr(1.0 - kg)) - (1.0 / sqr(x - kg))) /
    ((1.0 / sqr(1.0 - kg)) - (1.0 / sqr(kg)));
}



/*
 *  Actual shader!
 */
surface
strauss ( float Ka = 0.5,     
          t = 0.0,
          s = 0.7,
          m = 0.1;
        )
{
    float kj = 0.1;
    normal Nn = normalize(N); 
    float d = (1 - m*s);
    float h = 3.0 / (1.0 - s);
    float rd = (1.0 - s * s * s) * (1.0 - t); 
    float rn = 1.0 - t - rd;

    color Qa = rd * Cs;
    color Qd = (0.0, 0.0, 0.0);
    color Qs = (0.0, 0.0, 0.0);

    illuminance(P, Nn, PI/2.0)
    {
        normal Ln = normalize(L);
        float NdotL = max(Nn . Ln, 0.0);
        
        if(NdotL > 0.0)
        {
            Qd += Cl * NdotL * d * rd * Cs;

            normal V = normalize(I); 
            normal H = normalize(Ln - 2.0*NdotL*Nn);
            float HdotV = H . V;
            float NdotV = max(Nn . V, 0);
            float j = F(NdotL) * G(NdotL) * G(NdotV);
            float rj = min(1.0, rn + (rn + kj)*j);
            float rs = pow(HdotV, h) * rj;
            
            color CspecDiff = m * (1.0 - F(NdotL)) * (Cs - Cl);

            color Cspec = Cl + CspecDiff;

            //if(Cspec < (0.0, 0.0, 0.0))
            //    Cspec = Cl - CspecDiff; 
            
            Qs += Cl * rs;
        }
    }  
    
    Ci = Qs + Qd + Ka*Qa;
    Oi = Os;  
    Ci *= Oi;
}
